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Abstract 

This  paper  presents  the  results  of  a  parametric  study  conducted  with  a  previously  described  three-dimensional,  non-isothermal  model  of 
a  polymer  electrolyte  membrane  (PEM)  fuel  cell.  The  effect  of  various  operational  parameters  such  as  the  temperature  and  pressure  on  the 
fuel  cell  performance  was  investigated  in  detail.  It  was  found  that  in  order  to  obtain  physically  realistic  results  experimental  measurements 
of  various  modelling  parameters  were  needed.  The  results  show  good  qualitative  agreement  with  experimental  results  published  in  the 
literature.  In  addition,  geometrical  and  material  parameters  such  as  the  gas  diffusion  electrode  (GDE)  thickness  and  porosity  as  well  as  the 
ratio  between  the  channel  width  and  the  land  area  were  investigated.  The  contact  resistance  inside  the  cell  was  found  to  play  an  important 
role  for  the  evaluation  of  the  impact  of  such  parameters  on  the  fuel  cell  performance.  The  results  demonstrate  the  usefulness  of  this 
computational  model  as  a  design  and  optimization  tool. 

©  2003  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Fuel  cells  convert  the  chemical  energy  of  hydrogen  and 
oxygen  directly  into  electricity.  Their  high  efficiency  and 
low  emissions  have  made  them  a  prime  candidate  for  pow¬ 
ering  the  next  generation  of  electric  vehicles,  and  their  mod¬ 
ular  design  and  the  prospects  of  micro-scaling  them  have 
gained  the  attention  of  cellular  phone  and  laptop  manufac¬ 
turers.  Their  scalability  makes  them  prime  candidates  for  a 
variety  of  stationary  applications  including  distributed  resi¬ 
dential  power  generation.  The  basic  structure  and  operation 
principle  the  polymer  electrolyte  membrane  (PEM)  fuel  cell 
considered  here  are  illustrated  in  Fig.  1. 

The  polymer  electrolyte  consists  of  a  perfluorinated  poly¬ 
mer  backbone  with  sulfonic  acid  side  chains.  When  fully 
humidified,  this  material  becomes  an  excellent  protonic  con¬ 
ductor.  The  membrane  and  the  two  electrodes  (teflonated 
porous  carbon  paper  or  cloth  with  platinum  on  supported 
carbon)  are  assembled  into  a  sandwich  structure  to  form  a 
membrane-electrode  assembly  (ME  A).  The  ME  A  is  placed 
between  two  graphite  bipolar  plates  with  machined  groves 
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that  provide  flow  channels  for  distributing  the  fuel  (hydro¬ 
gen)  and  oxidant  (oxygen  from  air). 

The  hydrogen-rich  fuel  is  fed  to  the  anode,  where  the  hy¬ 
drogen  diffuses  through  the  porous  gas  diffusion  electrode 
(GDE).  At  the  catalyst  layer,  the  hydrogen  splits  into  hydro¬ 
gen  protons  and  electrons  according  to: 

2H2^4H++4e-  (1) 

Driven  by  an  electric  field,  the  H+  ions  migrate  through  the 
polymer  electrolyte  membrane.  The  oxygen  in  the  cathode 
gas  stream  diffuses  through  the  gas  diffusion  electrode  to¬ 
wards  the  catalyst  interface  where  it  combines  with  the  hy¬ 
drogen  protons  and  the  electrons  to  form  water  according  to: 

02  +  4H+  +  4e"  =►  2H20  (2) 

The  overall  reaction  is  exothermic  and  can  be  written  as: 

2H2  +  O2  =)>  2H2O  +  electricity  +  heat  (3) 

Several  coupled  fluid  flow,  heat  and  mass  transport  pro¬ 
cesses  occur  in  a  fuel  cell  in  conjunction  with  the  electro¬ 
chemical  reaction.  These  processes  have  a  significant  impact 
on  two  important  operational  issues:  (i)  thermal  and  wa¬ 
ter  management,  and  (ii)  mass  transport  limitations.  Water 
management  ensures  that  the  polymer  electrolyte  membrane 
remains  fully  hydrated  to  maintain  good  ionic  conductivity 
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Fig.  1.  Schematic  of  a  PEM  fuel  cell. 


and  performance.  Water  content  of  the  membrane  is  deter¬ 
mined  by  the  balance  between  water  production  and  three 
water  transport  processes:  electro-osmotic  drag  of  water, 
associated  with  proton  migration  through  the  membrane 
from  the  anode  to  the  cathode  side;  back  diffusion  from  the 
cathode;  and  diffusion  of  water  to/from  the  oxidant/fuel  gas 
streams.  Without  control,  an  imbalance  between  production 
and  removal  rates  of  water  can  occur.  This  results  in  either 
dehydration  of  the  membrane,  or  flooding  of  the  electrodes; 
both  phenomena  have  a  very  detrimental  effect  on  perfor¬ 
mance  and  fuel  cells  have  to  be  carefully  designed  to  avoid 
the  occurrence  of  either  phenomenon.  Thermal  management 
is  required  to  remove  the  heat  produced  by  the  electrochem¬ 
ical  reaction  in  order  to  prevent  drying  out  of  the  membrane 
and  excessive  thermal  stresses  that  may  result  in  rupture  of 
the  membrane.  The  small  temperature  differentials  between 
the  fuel  cell  stack  and  the  operating  environment  make 
thermal  management  a  challenging  problem  in  PEMFC’s. 

Because  of  the  highly  reactive  environment  and  compact 
nature  of  a  fuel  cell  it  is  not  possible  to  perform  detailed 
in  situ  measurements  during  operation.  Such  information 
has  been  sought  through  modelling  and  simulation  in  or¬ 
der  to  improve  understanding  of  water  and  species  trans¬ 
port,  optimize  thermal  management  and  shorten  the  design 
and  optimization  cycles.  Modelling  of  fuel  cells  is  chal¬ 
lenging,  because  the  processes  involve  multi-component, 
multi-phase,  and  multi-dimensional  flow,  heat  and  mass 
transfer  with  electrochemical  reactions  all  occurring  in 
irregular  geometries  including  porous  media.  Numerous 
authors  have  developed  fuel  cell  models  accounting  for  var¬ 
ious  physical  processes.  The  most  prominent  earlier  works 
stem  from  Bernardi  and  Verbrugge  [1,2]  and  Springer  et  al. 
[3],  who  developed  one-dimensional,  isothermal  models  of 
the  membrane-electrode  assembly.  Fuller  and  Newman  [4] 
published  a  quasi  two-dimensional  model  based  on  con¬ 
centration  theory.  The  work  by  Nguyen  and  White  [5],  Yi 
and  Nguyen  [6]  was  two-dimensional  in  nature,  but  the 


gas  diffusion  electrodes  were  omitted,  assuming  “ultrathin” 
electrodes.  The  importance  of  accounting  for  temperature 
gradients  in  fuel  cells  modelling  was  demonstrated  in  the 
work  of  Wohr  et  al.  [7]  and  Djilali  and  Fu  [8].  The  impor¬ 
tant  issue  of  water  flooding  was  addressed  by  Baschuk  and 
Fi  in  a  one-dimensional  model  [9]. 

Earlier  models  were  primarily  analytic  and  required  a 
number  of  simplifications  due  to  the  limitations  of  the  nu¬ 
merical  techniques.  More  recently,  a  general  trend  can  be 
observed  to  apply  the  methods  of  computational  fluid  dy¬ 
namics  to  fuel  cell  modelling.  Gurau  et  al.  [10]  published  a 
fully  two-dimensional  model  of  a  whole  fuel  cell,  i.e.  two 
gas  flow  channels  separated  by  the  membrane-electrode  as¬ 
sembly.  Wang  et  al.  [11,12]  have  developed  a  similar  model 
and  included  two-phase  flow. 

In  an  earlier  publication,  our  research  group  has  pre¬ 
sented  a  non-isothermal,  three-dimensional  model  of  a  com¬ 
plete  single  cell  [13].  Phase-change  and  multi-phase  flow, 
however,  could  not  be  addressed  with  that  model.  Recently, 
Dutta  et  al.  [14]  published  a  three-dimensional  computa¬ 
tional  model  based  on  the  commercial  software  package  Flu¬ 
ent ,  which  is  quite  similar  to  the  one  we  presented  earlier, 
except  that  it  accounts  for  a  partially  dehydrated  membrane 
using  an  empirical  approach. 

Overall,  it  can  be  said  that  whereas  the  focus  of  earlier 
modelling  efforts  was  often  one-dimensional  with  focus  on 
the  electrochemistry,  more  recent  publications  utilize  the 
methods  of  computational  fluid  dynamics  in  order  to  devise 
multi-dimensional  and  multi-phase  models. 

2.  Model  description 

A  detailed  description  of  the  model  has  recently  been  pub¬ 
lished  [13]  and  shall  not  be  repeated  here.  In  brief,  the  model 
is  based  on  the  commercial  software  package  CFX-4.3  (AEA 
Technology)  and  considers  single-phase,  multi-component 
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flow  inside  the  gas  flow  channels  and  the  porous  media.  Wa¬ 
ter  vapour  is  assumed  to  exist  at  its  saturation  pressure  in 
both  gas  streams.  The  physical  domain  consists  of  a  straight 
channel  of  5  cm  length  and  rectangular  cross-section.  Be¬ 
cause  of  the  symmetry  assumption  only  half  a  channel  and 
half  a  land  area  in  between  has  to  be  modelled.  For  com¬ 
putational  convenience  the  physical  domain  has  been  split 
up  into  four  domains,  allowing  for  conjugate  heat  transfer 
between  the  fluid-phase  and  the  solid  matrix  of  the  porous 
medium  in  one  domain,  liquid  water  flux  through  the  mem¬ 
brane  and  the  GDE  in  a  second  subdomain  and  the  electrical 
field  in  the  membrane  in  the  fourth  domain  with  appropriate 
boundary  conditions. 

Overall,  the  model  solves  for  the  multi-component  flow 
in  the  channels  and  the  porous  media  with  heat  transfer. 
In  addition,  several  phenomenological  equations  are  imple¬ 
mented  in  a  suite  of  user-subroutines,  which  are: 


conducted  with  this  model  in  order  to  assess  the  impact 
of  the  various  fuel  cell  operating  and  material  parame¬ 
ters  on  the  fuel  cell  performance.  In  the  following  section 
only  the  parameter  investigated  is  changed,  all  other  pa¬ 
rameters  are  at  the  base  case  conditions  as  outlined  in 
[13]. 


3.  Results  and  discussion 

3. 1.  Effect  of  temperature 

In  order  to  successfully  model  the  effect  of  the  tempera¬ 
ture  on  the  fuel  cell  performance,  a  basic  understanding  of 
its  direct  influence  on  various  model  parameters  is  required. 
The  properties  most  dependent  on  temperature  are  as  fol¬ 
lowing. 


•  the  generalized  Fick’ s  law  for  multi-component  diffusion 

[15]; 

•  the  Butler- Volmer  equation  for  reaction  kinetics  [16]; 

•  the  Nernst-Planck  equation  for  the  transport  of  protons 
through  the  membrane  [16];  and 

•  a  modified  version  of  the  Schlogl  equation  to  account  for 
the  flux  of  liquid  water  through  the  membrane  [17]. 


The  cell  voltage  is  calculated  at  the  end  of  each  computa¬ 
tion,  accounting  for  the  various  loss  mechanisms,  according 
to: 


?7act  ^contact  Omtm 


where  the  reversible  cell  potential  is  being  calculated 
as  function  of  the  temperature  and  pressure  using  the  Nernst 
equation,  and  the  activation  overpotential  ijnct  is  calculated 
employing  the  Tafel  equation  by  assuming  an  apparent  ex¬ 
change  current  density  that  depends  on  parameters  such  as 
the  operation  temperature,  pressure  and  the  catalyst  loading. 
In  contrast  to  [13]  in  the  current  paper  the  ohmic  loss  also 
includes  a  loss  due  to  contact  resistance  at  the  interface  be¬ 
tween  the  bipolar  plates  and  the  gas  diffusion  layers,  which 
have  so  far  been  neglected  in  virtually  all  modelling  studies. 
The  voltage  loss  due  to  contact  resistance  can  be  assumed 
of  ohmic  nature,  i.e. 


•  The  composition  of  the  incoming  gas  streams.  Assuming 
the  inlet  gases  are  fully  humidified,  the  partial  pressure  of 
water  vapour  entering  the  cell  depends  on  the  temperature 
only.  Thus,  the  molar  fraction  of  water  vapour  is  a  function 
of  the  inlet  pressure  and  temperature,  and  so  the  molar 
fraction  of  the  incoming  hydrogen  and  oxygen  depend  on 
the  temperature  and  pressure  as  well. 

•  The  exchange  current  density  z'o  of  the  oxygen  reduction 
reaction  (ORR)  increases  rapidly  with  temperature  due  to 
the  enhanced  reaction  kinetics.  Parthasarathy  et  al.  [18] 
conducted  experiments  in  order  to  determine  a  correlation 
between  the  cell  temperature  and  the  exchange  current 
density  of  the  oxygen  reduction  reaction. 

•  The  membrane  conductivity  k  increases,  because  a  higher 
temperature  leads  also  to  a  higher  diffusivity  of  the  hy¬ 
drogen  protons  in  the  electrolyte  membrane,  thereby  re¬ 
ducing  the  membrane  resistance  [1]. 

•  The  reference  potential  £°.  Although  the  Nernst  equation 
[16]  shows  a  decrease  in  the  reference  potential  with  an 
increasing  temperature,  experimental  results  indicate  an 
increase,  which  can  be  explained  with  a  higher  diffusivity 
of  the  hydrogen  with  increasing  temperature  [18]. 

•  The  gas-pair  diffusivities  Dij  in  the  Stefan-Maxwell  equa¬ 
tions.  An  increase  in  temperature  leads  to  an  increase  in 
the  binary  gas-pair  diffusivities  [19]. 


^contact  —  ^contact  (P/ 

where  rcontact  is  an  assumed  contact  resistance  in  (Q  cm2) 
that  varies  with  the  porosity  of  the  gas  diffusion  layer  and  the 
contact  area  between  the  bipolar  plates  and  the  gas  diffusion 
layer,  i.e.  the  ratio  between  the  channel  width  and  the  land 
area.  This  will  be  explained  in  detail,  later.  The  membrane 
loss  77mem  is  calculated  assuming  the  membrane  is  fully  hy¬ 
drated  at  all  times  so  that  the  protonic  conductivity  remains 
constant. 

The  basic  capabilities  and  limitations  of  the  current 
model  have  been  outlined  in  detail  in  an  earlier  publica¬ 
tion  [13].  The  current  paper  focuses  on  a  parametric  study 


In  order  to  determine  the  inlet  gas  composition  as  a  function 
of  temperature,  the  following  relation  between  the  tempera¬ 
ture  and  the  saturation  pressure  of  water  has  been  used  [3]: 

log10  psat  =  -2.1794  +  0.02953  x  §  -  9.1837E  -5xil2 
+  1 .4454E  -  7  x  tf3  (6) 


where  d  is  the  temperature  (°C).  The  molar  fraction  of  water 
vapour  in  the  incoming  gas  stream  is  simply  the  ratio  of  the 
saturation  pressure  and  the  total  pressure: 


7*  sat 

-*H20,in  =  - 

Pin 
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Fig.  2.  Molar  inlet  composition  of  the  cathode  side  gas  stream  as  function  of  temperature. 


Since  the  ratio  of  nitrogen  and  oxygen  in  dry  air  is  known 
to  be  79:21,  the  inlet  oxygen  fraction  can  be  found  via: 

_  1  -  *H2Q,in 

°2’in  i  +  (79/21)  } 

The  resulting  inlet  gas  composition  for  different  pressures  is 
shown  in  Fig.  2.  Clearly,  at  an  operating  pressure  of  1  atm  the 
effect  of  the  temperature  on  the  inlet  composition  is  much 
stronger  than  at  elevated  pressures.  At  80  °C  for  atmospheric 
pressure,  almost  50%  (molar)  of  the  incoming  cathode  side 
gas  stream  consists  of  water  vapour  and  only  around  1 1  % 
is  oxygen. 

In  order  to  find  a  correlation  between  the  reference  ex¬ 
change  current  density  io  for  the  oxygen  reduction  reaction 
at  the  cathode  side  and  the  temperature,  experimental  re¬ 
sults  obtained  by  Parthasarathy  et  al.  [18]  were  used.  The 
catalyst  loading  in  those  experiments  did  not  correspond  to 
the  loading  in  our  simulated  cell.  A  method  for  estimating 
the  variation  of  the  exchange  current  density  with  tempera¬ 
ture  was  devised  (see  [20])  and  is  used  here.  Table  1  shows 
the  resulting  exchange  current  densities  used  in  our  model 
for  a  range  of  temperatures.  The  exchange  current  density 
for  a  given  catalyst  loading  varies  by  an  order  of  magnitude 
from  323  to  353  K,  and  this  is  consistent  with  the  results  of 
Parthasarathy  et  al.  [18].  Note  that  the  base  case  temperature 
in  this  study  was  353  K. 

The  conductivity  of  the  electrolyte  membrane  is  also  a 
strong  function  of  temperature.  A  theoretical  expression  de¬ 
rived  by  Bernardi  and  Verbrugge  [2]  shows  that  is  depends 


linearly  on  the  diffusivity  of  the  hydrogen  proton  inside  the 
membrane  Z)H+ ,  which,  in  turn,  was  calculated  according  to 
[21]: 


DH+fi(T) 

- =  const. 

T 


where  /x(7)  was  taken  by  Bernardi  and  Verbrugge  to  be  the 
viscosity  of  water.  The  resulting  diffusivities  for  various  tem¬ 
peratures  are  listed  in  the  second  column  of  Table  2,  using 
the  value  of  1.4  x  10-5  cm2/s  at  22  °C  as  measured  by  Ver¬ 
brugge  and  Hill  [22]  as  a  benchmark.  However,  membrane 
conductivities  computed  using  this  approach  do  not  match 
measured  values.  For  the  present  simulation  the  membrane 
conductivity  was  taken  to  be  6.8  x  1CT2  S/cm  at  80  °C  [3], 
and  a  linear  scaling  with  protonic  diffusivity  was  used  to 
estimate  values  at  other  temperatures. 

The  reversible  cell  potential  at  various  temperatures 
was  computed  using  the  Nernst  equation: 

RT 

E°  =  1.23  -  0.9  x  10“3(T  -  298)  +  2.3  —  pl,po2  (10) 

4  F  2 


Finally,  the  binary  diffusivities  of  the  gas-phase  pairs  had  to 
be  scaled  with  temperature.  Different  theoretical  predictions 
of  the  binary  diffusivity  values  can  be  found  in  the  literature, 
and  the  one  taken  here  was  [19]: 


Dij{T)  =  Ay  (7o) 


T 

7b 


1.75 


(ID 


The  polarization  curves  obtained  for  various  cell  temper¬ 
atures  and  using  the  above  parameter  values  are  shown  in 


Table  1 

Apparent  exchange  current  density  of  the  ORR  as  a  function  of  temper¬ 
ature 


Table  2 

Proton  diffusivity  and  membrane  conductivity  as  function  of  temperature 


T  (K) 

io  (A/cm2) 

T( K) 

Dh+  (cm2/s) 

k  (S/cm) 

323 

3.3  x  10"8 

353 

4.54  x  10“5 

0.068 

333 

7.9  x  IO"7 

343 

3.87  x  10“5 

0.058 

343 

1.7  x  IO"7 

333 

3.26  x  10“5 

0.049 

353 

4.4  x  10“7 

323 

2.75  x  10“5 

0.041 
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(a)  Current  Density  [A/cm2]  (b)  Current  Density  [A/cm2] 

Fig.  3.  Predicted  polarization  curves  (a)  and  power  density  curves  (b)  for  different  cell  temperatures.  All  other  parameters  are  at  base  case  conditions. 


Fig.  3.  The  change  in  the  initial  drop  due  to  the  lower  ex¬ 
change  current  density  is  relatively  small  compared  to  the 
drop-off  in  the  linear  region,  where  ohmic  losses  are  pre¬ 
dominant.  The  limiting  current  density  is  very  similar  in  all 
cases.  This  is  due  to  the  fact  that  the  base  case  pressure  is 
high  (5  atm),  and  hence  the  inlet  gas  composition  changes 
little  with  temperature  (cf.  Fig.  2).  Changes  in  operating 
pressure  on  the  other  hand  have  a  large  impact  on  the  inlet 
composition  and,  hence,  on  the  limiting  current  density,  as 
will  be  shown  below.  We  note  in  Fig.  3  that  the  maximum 
power  density  shifts  towards  higher  current  density  with  an 
increasing  temperature  as  a  result  of  reduced  ohmic  losses 
[23]. 

A  systematic  comparison  with  experimental  results  is  dif¬ 
ficult,  since  the  exact  experimental  conditions  are  frequently 
incompletely  documented.  As  the  current  model  has  histor¬ 
ically  been  developed  as  a  three-dimensional  version  of  the 
model  by  Bernardi  and  Verbrugge  [2],  who  in  turn  matched 
their  modelling  data  to  experiments  conducted  by  Ticianelli 
et  al.  [24],  this  shall  be  the  first  comparison  made.  Fig.  4 


Current  Density  [A/cm2] 

Fig.  4.  Polarization  curves  for  different  temperatures  obtained  by  Ticianelli 
et  al. 


shows  experimental  results  obtained  for  two  different  fuel 
cells  and  at  different  temperatures.  Qualitatively,  the  results 
agree  with  our  numerical  data. 

Kim  et  al.  [25]  extended  the  curve-fitting  approach  that 
was  suggested  by  Ticianelli  et  al.  in  order  to  include  the 
drop-off  due  to  mass  transport  limitations  at  high  current 
densities.  The  corrected  equation  is  [25]: 

E  =  Er  —  blog- - Rfi  —  m  exp(m)  (12) 

where  Ev  is  the  reversible  cell  potential  for  the  given  condi¬ 
tions  and  the  three  following  terms  describe  the  various  loss 
mechanisms.  The  first  term  can  be  recognized  as  the  Tafel 
equation  that  describes  the  activation  overpotential,  which 
is  predominant  at  low  current  densities.  The  second  term  Ri 
describes  a  linear  drop-off,  which  is  predominant  in  the  in¬ 
termediate  current  density  region,  where  Ri  is  the  internal 
resistance  caused  by  membrane  and  contact  losses.  The  last 
term  becomes  predominant  in  the  high  current  density  re¬ 
gion,  and  is  used  to  match  the  drop-off  towards  the  limiting 
current  density.  A  physical  interpretation  for  the  parame¬ 
ters  m  and  n  was  not  given,  but  Bevers  et  al.  [26]  found  in 
their  one-dimensional  modelling  study  that  m  correlates  to 
the  electrolyte  conductivity  and  n  to  the  porosity  of  the  gas 
diffusion  layer.  Following  up  on  this  we  can  speculate  now 
that  both  m  and  n  relate  to  water  management  issues:  a  par¬ 
tially  dehydrated  electrolyte  membrane  leads  to  a  decrease 
in  conductivity,  which  can  be  represented  by  m,  whereas  an 
excess  in  liquid  water  leads  to  a  reduction  in  porosity  and 
hence  to  an  early  onset  of  mass  transport  limitations,  which 
can  be  captured  by  the  parameter  n.  Recently,  Natarajan  and 
Nguyen  [27]  have  shown,  how  strongly  the  cell  potential  de¬ 
pends  on  the  saturation  level  of  the  air.  Note  also  that  both 
these  effects  are  not  accounted  for  by  the  first  two  loss  terms 
in  Eq.  (12). 

Fig.  5  shows  the  effect  of  temperature  on  the  fuel  cell  per¬ 
formance,  as  was  measured  by  Kim  et  al.  [25].  Note  that  in 
this  case  the  stoichiometric  flow  ratio  was  only  1.5,  which 
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Fig.  5.  Polarization  curves  for  different  temperatures  and  pressures  repro¬ 
duced  from  Kim  et  al. 

led  to  relatively  low  limiting  current  densities.  The  impor¬ 
tance  of  the  stoichiometric  flow  ratio  is  discussed  in  detail  in 
a  later  section.  The  qualitative  agreement  between  the  exper¬ 
imental  data  and  numerical  predictions  is  good,  except  for 
the  high  current  density  region.  This  is  due  to  the  fact  that 
the  calculated  cell  voltage  in  our  model  assumes  a  fully  hy¬ 
drated  membrane  at  all  times  and  a  constant  porosity  of  the 
gas  diffusion  layers,  neglecting  the  onset  of  pore-plugging 
by  liquid  water. 

3.2.  Effect  of  pressure 

Similar  to  the  temperature,  the  operating  pressure  affects 
numerous  parameter  that  are  important  for  the  fuel  cell  op¬ 
eration,  among  which  are  as  follows. 

•  The  inlet  gas  compositions.  Because  the  saturation  pres¬ 
sure  of  water  vapour  depends  only  on  the  temperature,  a 
change  in  the  operating  pressure  leads  to  a  change  in  the 
inlet  gas  compositions,  assuming  the  inlet  gases  are  fully 
humidified. 

•  The  exchange  current  density  io  of  the  oxygen  reduction 
reaction.  The  dependence  of  the  cathodic  exchange  cur¬ 
rent  density  on  the  oxygen  pressure  was  also  investigated 
experimentally  by  Parthasarathy  et  al.  [28]. 

•  The  reference  potential  Ev.  According  to  the  Nernst  equa¬ 
tion,  an  increased  pressure  leads  to  an  increase  in  the  equi¬ 
librium  potential. 

•  The  gas-pair  diffusivities  Dij  in  the  Stefan-Max  well  equa¬ 
tions.  It  is  well  known  that  the  product  of  pressure  and 
the  binary  diffusivity  is  constant  [19].  Hence,  a  doubling 
of  the  pressure  will  cut  the  binary  diffusivity  in  half. 

Since  the  saturation  pressure  for  water  is  only  a  function 
of  temperature,  it  remains  constant  for  a  variation  of  the 
inlet  pressure,  and  the  molar  fraction  of  water  vapour  in 
the  incoming  cathode  gas  stream  is  given  by  Eqs.  (6)  and 
(7).  The  molar  oxygen  fraction  results  then  out  of  Eq.  (8). 


log  [P02/atm] 


Fig.  6.  The  dependence  of  the  exchange  current  density  of  the  oxygen 
reduction  reaction  on  the  oxygen  pressure. 

It  was  already  noted  in  Fig.  2  that  the  change  in  the  inlet 
gas  composition  is  particularly  strong  in  the  range  from  1  to 
3  atm.  Above  3  atm  the  composition  changes  only  slightly 
with  the  pressure. 

In  order  to  understand  the  dependence  of  the  ex¬ 
change  current  density  io  on  the  partial  oxygen  pressure, 
Parthasarathy  et  al.  [28]  conducted  experiments  at  a  tem¬ 
perature  of  50  °C.  The  results  are  summarized  in  Fig.  6. 

A  linear  relationship  was  found  between  the  logarithm 
of  the  exchange  current  density  io  and  the  logarithm  of  the 
oxygen  partial  pressure,  according  to: 

io  =  1.27  x  10-8  exp(2.06/?o2)  (13) 

Because  the  temperature  was  different  from  our  case  and 
the  roughness  factor  or  catalyst  loading  was  not  reported, 
the  values  obtained  with  Eq.  (12)  had  to  be  further  scaled 
using  the  base  case  value  of  io  =  4.4  x  10-7  A/cm2  as  a 
reference  value.  The  detailed  procedure  is  outlined  in  [20]. 

Table  3  summarizes  the  exchange  current  density  used 
for  our  simulations.  The  third  column  lists  the  exchange 
current  density  as  was  calculated  using  Eq.  (12)  and  the 
partial  oxygen  pressure  as  listed  in  the  second  column,  and 
the  last  column  contains  the  scaled  values  for  the  different 
temperature  and  catalyst  loading.  It  can  be  seen  that  the 
pressure  has  a  much  weaker  effect  on  the  exchange  current 
density  than  the  temperature. 

Although  it  is  known  that  the  reference  exchange  potential 
at  the  cathode  side  is  a  mixed  potential  due  to  competing 

Table  3 


Exchange  current  density  of  the  ORR  as  a  function  of  pressure 


p  (atm) 

Po2  (atm) 

io 

T  =  50°C 

T  =  80 

°C 

1.0 

0.1251 

1.64  x  10"8 

0.78  x 

10-7 

1.5 

0.1483 

1.72  x  10“8 

0.82  x 

10-7 

3.0 

0.5451 

3.90  x  10“8 

1.85  x 

10“7 

5.0 

0.9650 

9.27  x  10"8 

4.40  x 

10-7 
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Fig.  7.  Average  molar  oxygen  fraction  at  the  cathode  side  catalyst  layer  ^g-  9.  Experimentally  observed  polarization  curves  for  various  cathode 

as  function  of  the  operating  pressure.  side  Pressures  at  a  temperature  of  70  °C  and  a  stoichiometric  flow  ratio 

of  1.5. 


reactions  [28],  the  adjustment  of  the  reference  potential  EP 
was  done  according  to  the  corrected  Nernst  equation  in  our 
case,  and  the  diffusion  coefficients  for  the  Stefan-Maxwell 
equations  were  adjusted  automatically  in  our  model. 

With  the  functional  variations  of  the  transport  parameters 
determined,  computations  were  performed  for  various  op¬ 
erating  pressures.  Fig.  7  shows  the  average  molar  oxygen 
fraction  at  the  cathode  side  catalyst  layer,  which  eventually 
determines  the  limiting  current  density.  The  stoichiometric 
flow  ratio  is  maintained  constant  at  3.0  in  all  these  cases. 
The  higher  oxygen  fraction  at  the  cathode  side  inlet  even¬ 
tually  leads  to  a  higher  maximum  current.  The  increase  is 
significant,  when  the  pressure  increases  from  atmospheric 
pressure  to  3  atm,  and  this  is  consistant  with  Fig.  2.  Further 
increase  in  pressure  from  3.0  to  5.0  atm  does  not  lead  to 
a  significant  improvement  in  terms  of  the  limiting  current 
density.  It  should  be  emphasized  again  that  this  presumes 
the  incoming  gases  are  fully  humidified. 

The  polarization  curves  shown  in  Fig.  8  reveal  a  signif¬ 
icant  change  in  the  initial  drop-off,  when  the  pressure  is 
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Fig.  8.  Calculated  polarization  curves  for  various  operating  pressures. 


changed.  This  can  be  attributed  to  the  change  in  the  equi¬ 
librium  potential  that  goes  along  with  a  decrease  in  the  re¬ 
actant  pressure  (Nernst  equation).  To  a  much  lesser  extend, 
the  decrease  in  the  exchange  current  density  with  decreasing 
pressure  also  contributes  to  this  effect. 

Again,  a  detailed  comparison  with  experimental  results 
from  the  literature  can  only  be  made  on  a  qualitative  basis, 
since  the  exact  conditions  of  the  various  experiments  are 
not  reported.  In  Fig.  9,  experimentally  obtained  polarization 
curves  by  Kim  et  al.  [25]  are  reproduced.  The  experiments 
were  conducted  with  pure  hydrogen  at  the  anode  side  and 
air  at  the  cathode  side.  Although  the  exact  details  of  the 
experiments  such  as  the  cell  geometry  are  not  known,  the 
two  main  effects  that  the  cathode  side  pressure  has  on  the 
fuel  cell  performance  can  be  observed:  the  increase  of  the 
limiting  current  density  with  an  increase  in  pressure  and  an 
overall  better  cell  performance,  which  can  be  attributed  to 
an  increase  in  the  equilibrium  potential  as  well  as  reduced 
activation  losses.  It  is  interesting  to  note  that  the  measured 
polarization  curves  for  3.0  and  5.0  atm  are  relatively  close, 
which  agrees  well  with  the  modelling  results  shown  above. 

In  general  it  is  difficult  to  compare  the  results  obtained 
with  the  current  model  with  experimental  results  taken  from 
the  literature,  since  various  parameters  that  are  not  given  in 
the  literature  influence  the  fuel  cell  performance.  The  grad¬ 
ual  decrease  in  performance  with  current  density  that  has 
been  observed  by  Kim  et  al.  [25]  can  currently  not  be  cap¬ 
tured  with  our  model  and  we  have  reason  to  believe  that  it  is 
associated  with  water  management  issues  such  as  a  partial 
dehydration  of  the  membrane.  Qualitative  agreement,  how¬ 
ever,  is  very  good  and  the  principal  physical  benefits  of  oper¬ 
ating  a  fuel  cell  at  an  elevated  pressure  have  been  confirmed. 

3.3.  Effect  of  stoichiometric  flow  ratio 

The  stoichiometric  flow  ratio  is  an  important  parameter, 
as  it  is  the  inverse  of  the  fuel  utilization  and  hence  directly 
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Fig.  10.  Average  molar  oxygen  fraction  at  the  cathode  side  catalyst  layer 
as  function  of  current  density  for  different  values  of  the  stoichiometric 
flow  ratio. 

affects  the  theoretical  fuel  cell  efficiency  [29].  In  practice, 
the  stoichiometric  flow  ratio  has  an  impact  on  the  limiting 
current  density,  but,  more  importantly,  on  water  manage¬ 
ment  issues  [29],  which  can  not  be  resolved  with  the  current 
version  of  our  model.  The  effect  of  the  stoichiometric  flow 
ratio  on  the  limiting  current  density  is  shown  in  Fig.  10. 
The  limiting  current  density  is  reached,  when  the  oxygen 
concentration  at  the  catalyst  layer  becomes  zero.  Because 
the  stoichiometric  flow  ratio  is  constant  even  for  extremely 
small  current  densities,  the  average  molar  oxygen  fraction  at 


the  catalyst  layer  approaches  the  inlet  molar  oxygen  fraction 
of  around  18%  only  for  very  high  stoichiometric  flow  ratios 
(because  the  currents  can  be  extremely  small,  stoichiomet¬ 
ric  flow  ratios  are  often  in  the  several  hundreds  in  this  re¬ 
gion).  Apart  from  that  it  can  be  observed  that  the  lines  are 
equidistant  to  one  another  and  an  increase  in  the  stoichio¬ 
metric  flow  ratio  leads  to  an  increase  in  the  limiting  current 
density  region.  However,  beyond  a  value  of  3.0  the  gain  in 
current  density  is  relatively  small,  compared  to  the  gain  in 
the  region  below  a  value  of  2.0.  This  graph  also  suggests 
that  there  is  an  inherent  limitation  to  the  attainable  current 
density,  as  the  slope  does  not  change  for  a  change  in  the 
stoichiometric  flow  ratio.  We  will  see  below  that  this  slope 
depends  on  various  parameters  such  as  the  GDE  porosity 
and  channel  spacing. 

A  higher  stoichiometric  flow  ratio  also  results  in  a  more 
even  distribution  of  the  local  current  density,  as  is  shown  in 
Fig.  1 1.  It  was  observed  before  that  under  the  conditions  in¬ 
vestigated  the  amount  of  current  generated  under  the  chan¬ 
nel  area  increases  almost  linearly  with  the  current  density, 
until  the  limiting  current  is  reached  [13].  This  also  leads  to 
the  fact  that  for  a  lower  stoichiometric  flow  ratio  at  a  con¬ 
stant  current  density,  there  is  a  much  stronger  distribution  of 
current  inside  the  cell,  the  maximum  local  current  density 
being  at  the  inlet  under  the  channel  area. 

3.4.  Effect  of  GDE  porosity 

The  porosity  of  the  gas  diffusion  layer  has  two  competing 
effects  on  the  fuel  cell  performance:  as  the  porous  region 
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provides  the  space  for  the  reactants  to  diffuse  towards  the 
catalyst  region,  an  increase  in  the  porosity  means  that  the 
onset  of  mass  transport  limitations  occurs  at  higher  current 
densities,  i.e.  it  leads  to  higher  limiting  currents.  The  ad¬ 
verse  effect  of  a  high  porosity  is  an  expected  increase  in  the 
contact  resistance.  Contact  resistance  occurs  at  all  interfaces 
inside  the  fuel  cell,  the  most  important  one  being  the  inter¬ 
face  between  the  bipolar  plates  and  the  gas  diffusion  lay¬ 
ers,  and  its  magnitude  depends  on  various  parameters  such 
as  the  surface  material  and  treatment  and  the  applied  stack 
pressure.  As  the  electrons  travel  through  the  solid  matrix  of 
the  GDE,  it  can  be  assumed  that  the  contact  losses  increase 
linearly  with  an  increased  porosity.  The  extend  of  contact 
losses  has  been  measured  by  Barbir  et  al.  [30],  who  found 
that  for  their  “standard”  fuel  cell,  the  contact  resistance  can 
be  as  high  as  150  m^  cm,  i.e.  the  voltage  loss  due  to  con¬ 
tact  resistance  at  a  current  density  of  1.0  A/cm2  would  be 
as  high  as  0.15  V.  In  order  to  include  the  loss  due  to  contact 
resistance  in  our  model,  we  assumed  a  resistance  at  base 
case  conditions  (s  =  0.4)  and  scaled  this  linearly  with  the 
porosity  of  the  GDE. 

Fig.  12  shows  the  average  molar  oxygen  fraction  at  the 
cathode  side  catalyst  layer  as  function  of  current  density.  It 
can  be  seen  that  the  porosity  has  a  significant  effect  on  the 
limiting  current  density.  An  increase  in  the  porosity  of  the 
GDE  from  £  =  0.3  to  0.5  more  than  triples  the  value  for  the 
attainable  current  density  from  around  0.75  A/cm2  to  just 
below  2.4  A/cm2.  This  is  because  the  transport  inside  the 
gas  diffusion  electrodes  is  primarily  diffusive,  and  the  binary 
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Fig.  12.  Average  molar  oxygen  fraction  at  the  catalyst  interface  as  function 
of  current  density  for  three  different  GDE  porosities. 
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gas  diffusivities  are  scaled  by  the  so-called  Bruggemann 
correction  [1],  where  the  porosity  affects  the  diffusivity  by 
a  power  of  1.5. 

In  addition,  a  higher  porosity  evens  out  the  local  current 
density  distribution,  as  shown  in  Fig.  13.  For  a  lower  value 
of  the  porosity  a  much  higher  fraction  of  the  total  current 
is  generated  under  the  channel  area.  This  can  lead  to  local 
hot-spots  inside  the  membrane,  which  can  lead  to  a  further 
drying  out  of  the  membrane,  thus  increasing  the  electric 
resistance,  which  in  turn  leads  to  more  heat  generation  and 
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Fig.  13.  Local  current  density  distribution  at  the  cathode  side  catalyst  layer  for  a  GDE  porosity  of  0.4  (a)  and  0.6  (b).  The  average  current  density  is 
1.0  A/cm2. 
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(a)  Current  Density  [A/cm2] 


Fig.  14.  Power  density  curves  for  three  different  GDE  porosities  and  an  assumed  contact  resistance  at  base  case  (s  =  0.4)  of  50m£2cm2  (a)  and 
150m£2cm2  (b). 


can  lead  to  a  failure  of  the  membrane.  Thus,  it  is  important 
to  keep  the  current  density  relatively  even  throughout  the 
cell. 

Fig.  14  shows  the  calculated  power  density  curves  for 
different  porosities  and  two  assumed  values  for  the  contact 
resistance  at  base  case  conditions.  It  can  be  seen  that  for  an 
assumed  contact  resistance  of  50  m£2  cm2  all  three  power 
density  curves  remain  fairly  close. 

As  was  seen  before,  the  curve  is  cut  off  for  a  porosity 
of  s  =  0.3  because  of  the  early  onset  of  mass  transport 
limitations.  For  a  porosity  of  £  =  0.5  the  maximum  attain¬ 
able  is  not  limited  by  the  mass  transport  anymore,  because 
this  would  occur  only  at  a  current  density  of  2.4  A/cm2. 
The  ohmic  loss  becomes  so  high  that  the  cell  voltage  and 
hence  the  power  density  become  zero  at  a  current  density 
of  1.8A/cm2.  For  a  contact  resistance  of  150  m£2  cm2  the 
ohmic  losses  at  the  highly  porous  GDE  become  so  high  that 
the  maximum  power  density  is  achieved  by  the  cell  with 
the  lowest  porosity.  However,  due  to  the  early  onset  of  mass 
transport  limitations,  a  porosity  of  0.4  appears  to  be  the  op¬ 
timum  under  the  current  conditions. 

3.5.  Effect  of  the  channel  width 

For  the  width  of  the  fuel  cell  channel  the  same  argu¬ 
ments  hold  as  for  the  porosity,  as  the  predominant  param¬ 
eters  affected  are  the  limiting  current  density  in  the  form 
of  mass  transport  limitations  and  on  the  other  hand  con¬ 
tact  resistance.  In  addition,  the  pressure  drop  inside  the  cell 
will  depend  on  the  channel  width.  In  our  computations,  we 
have  left  the  overall  pitch  between  channels  constant,  i.e. 
we  changed  the  ratio  between  the  channel  width  and  the 
shoulder  width,  which  in  the  base  case  is  one,  as  both  the 
channel  and  the  land  area  are  both  1  mm  wide.  Fig.  15 
shows  that  the  effect  of  the  channel  width  on  the  limiting 
current  density  is  not  quite  as  strong  as  the  effect  of  the 
porosity. 


The  local  current  density  distribution  for  the  cases  with 
a  channel  and  land  area  width  of  0.8  mm/ 1.2  mm  and 
1.2  mm/0.8  mm,  respectively,  is  shown  in  Fig.  16. 

For  the  narrow  channel  the  local  current  density  can  ex¬ 
ceed  more  than  twice  the  value  as  the  average  current  den¬ 
sity  with  a  sharp  drop-off  under  the  land  area,  where  the 
local  current  density  is  below  0.2  A/cm2.  The  wider  channel 
makes  for  a  much  more  evenly  distributed  current  through¬ 
out  the  cell. 

The  drawback  of  a  wider  channel  is  the  reduction  in  con¬ 
tact  area  between  the  GDE  and  the  bipolar  plates,  increas¬ 
ing  the  contact  resistance.  Fig.  17  shows  the  power  density 
curves  for  two  assumed  values  for  the  contact  resistance  at 
base  case  (Ch/L:  1  mm/1  mm).  The  advantage  in  terms  of 
power  density  becomes  pronounced  at  contact  resistance  of 
higher  than  50m£2cm2.  For  an  assumed  contact  resistance 
of  150  m^  cm2  at  the  base  case,  the  limiting  current  density 
is  the  same  for  the  cases  with  a  1  mm  wide  channel  and  a 
channel  width  of  1.2  mm,  because  at  the  widest  channel  the 


Fig.  15.  Average  molar  oxygen  fraction  at  the  cathode  side  catalyst  layer 
for  three  different  values  of  the  channel  and  land  area  width. 
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Fig.  16.  Local  current  density  distribution  at  the  cathode  side  catalyst  layer  for  a  channel  width  of  0.8  mm  (a)  and  1.2  mm  (b).  The  average  current 
density  is  1.0  A/cm2  in  both  cases. 


limiting  current  density  is  again  determined  by  the  ohmic 
losses,  not  by  mass  transport  limitations. 

3.6.  Effect  of  GDE  thickness 

The  effect  of  the  GDE  thickness  on  the  fuel  cell  perfor¬ 
mance  is  again  mostly  on  the  mass  transport,  as  the  ohmic 
losses  of  the  electrons  inside  the  GDE  can  be  neglected  due 
to  the  high  conductivity  of  the  carbon  fiber  paper.  Fig.  18 
shows  the  impact  of  the  GDE  thickness  on  the  attainable 
current  density.  It  is  also  interesting  to  note  that  at  low  cur¬ 


(a)  Current  Density  [A/cm2] 


rent  densities  an  increased  GDE  thickness  leads  to  a  higher 
average  oxygen  fraction  than  for  a  thinner  GDE.  The  reason 
for  this  is  that  a  thicker  GDE  gives  the  oxygen  more  room 
to  spread  out  in  the  lateral  (z-)  direction  under  the  land  area. 
The  current  density,  where  these  lines  cross,  depends  on  the 
ratio  between  the  channel  width  and  the  land  area,  and  it  is 
worth  exploring,  how  strongly.  The  right  hand  side  of  Fig.  18 
shows  the  same  plot  for  a  channel  width  of  0.8  mm  and  a 
land  area  width  of  1.2  mm.  In  this  case,  the  lines  cross  at  a 
current  density  of  0.4  A/cm2,  which  means  that  this  effect 
is  not  of  practical  importance. 


(b)  Current  Density  [A/cm2] 


Fig.  17.  Power  density  curves  for  different  channel  widths  and  two  assumed  values  for  the  contact  resistance  at  base  case. 
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Fig.  18.  Average  molar  oxygen  fraction  at  the  cathode  side  catalyst  layer  for  different  GDE  thicknesses  as  function  of  current  density.  Part  (a)  is  for  a 
channel  width  of  1mm  and  (b)  for  0.8  mm. 


4.  Conclusions 

A  parametric  study  using  a  three-dimensional  model  of 
a  PEM  fuel  cell  has  been  performed.  The  study  quantifies 
the  impact  of  operating  and  geometric  parameters  on  per¬ 
formance.  Experiments  with  respect  to  the  effect  of  tem¬ 
perature  and  pressure  are  well  documented  in  the  literature, 
and  the  results  obtained  with  our  model  reproduce  all  the 
trends  and  highlight  the  importance  of  properly  accounting 
for  the  effect  of  operational  parameters  on  transport  and  ki¬ 
netic  properties. 

In  addition,  material  properties  such  as  the  GDE  thick¬ 
ness  and  porosity  as  well  as  the  channel  width  compared  to 
the  land  area  have  been  investigated.  It  was  found  that  the 
porosity  of  the  gas  diffusion  layer  has  a  strong  effect  on  the 
limiting  current  density.  In  order  to  properly  assess  the  im¬ 
pact  of  porosity  and  channel  width  on  performance,  it  was 
necessary  to  estimate  the  extend  of  contact  resistance  inside 
the  fuel  cell. 

Overall,  this  study  demonstrates  the  use  of  a  comprehen¬ 
sive  three-dimensional  single-phase  model  as  a  design  tool. 
Further  capability  enhancement  requires  the  modelling  of 
two-phase  transport  and  phase  change  in  the  gas  diffusion 
electrodes.  Work  is  in  progress  that  addresses  these  issues. 
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